Radiation detection with non-parametric decompounding of pulse pile-up

ABSTRACT

A method of determining a spectrum of energies of individual quanta of radiation received in a radiation detector is disclosed. Spectrum sensitive statistics are computed from a time series of digital observations from the radiation detector, defining a mapping from a density of amplitudes of the pulses to the spectrum sensitive statistics. The spectrum is determined by estimating the density of amplitudes of the pulses by applying an inversion of the mapping to the spectrum sensitive statistics. The statistics may be based on a first set of nonoverlapping time intervals of constant length L at least as long as a duration of the pulses without regard to entirety of clusters of the pulses; and a second set of nonoverlapping time intervals of constant length L1 less than L also without regard to entirety of clusters of the pulses. A method of estimating count rate is also disclosed.

FIELD OF THE INVENTION

A novel method is provided for estimating the energy distribution of quanta of radiation such as photons incident upon a detector in spectroscopic systems such as in X-ray or gamma-ray spectroscopy. The method is particularly useful for count-rate regimes where pulse pile-up is an issue. A key step in the derivation of the estimator of an embodiment of the invention is the novel reformulation of the problem as a decompounding problem of a compound Poisson process. The method can be applied to any form of radiation detector detecting quanta or other corpuscles of radiation, such as x-rays, gamma rays or other photons, neutrons, atoms, molecules, or seismic pulses. Applications of spectroscopy from such detectors are well known. Such applications are described widely in the prior art, including in international patent applications PCT/AU2005/001423, PCT/AU2009/000393, PCT/AU2009/000394, PCTAU2009/000395, PCT/AU2009/001648, PCT/AU2012/000678, PCT/AU2014/050420, PCT/AU2015/050752, PCT/AU2017/050514, and PCT/AU2017/050512, each of which is incorporated herein in its entirety for the purpose of describing potential applications of the current invention and any other background material needed to understand the current invention.

1 BACKGROUND

X-RAY and gamma-ray spectroscopy underpin a wide range of scientific, industrial and commercial processes. One goal of spectroscopy is to estimate the energy distribution of photons incident upon a detector. From a signal processing perspective, the challenge is to convert the stream of pulses output by a detector into a histogram of the area under each pulse. Pulses are generated according to a Poisson distribution whose rate corresponds to the intensity of X-rays or gamma-rays used to illuminate the sample. Increasing the intensity results in more pulses per second on average and hence a shorter time before an accurate histogram is obtained. In applications such as baggage scanning at airports, this translates directly into greater throughput. Pulse pile-up occurs when two or more pulses overlap in the time domain. As the count rate (the average number of pulses per second) increases so does the incidence of pulse pile-up. This increases the difficulty in determining the number of pulses present and the area under each pulse. In the limit, the problem is ill-conditioned: if two pulses start at essentially the same time, their superposition is indistinguishable from a single pulse. The response of an X-ray or gamma-ray detector to incident photons can be modelled as the superposition of convolutions of pulse shape functions Φ_(j)(t) with Dirac-delta impulses,

$\begin{matrix} {{r(t)} = {\sum\limits_{j = {- \infty}}^{\infty}{a_{j}{\delta\left( {t - \tau_{j}} \right)}{{{\bigstar\Phi}_{j}(t)}.}}}} & (1) \end{matrix}$

The arrival times . . . , τ⁻¹, τ₀, τ₁, . . . are unknown and form a Poisson process. Each photon arrival is modelled as a Dirac-delta at time τ_(j), and with amplitude a_(j) that is proportional to the photon energy and induces a detector pulse shape response Φ_(j). The amplitudes a_(j) are realizations of identically distributed random variables A_(j) whose common probability density function ƒ_(A)(x) is unknown. The pulse shape function Φ_(j) is determined by the geometry of the detector and the photon interaction. In some systems the variation in pulse shape is minimal and may be ignored, while other systems (e.g., HPGe) individual pulse shapes may differ significantly from one another ill. It is assumed all pulse shapes are causal i.e., Φ_(j)(t)=0 for t<0, uni-modal, of finite energy, and decay exponentially toward zero as t→∞. The integral of the pulse shape functions are normalized to unity, i.e., ∫_(−∞) ^(∞)Φ_(j)(t) dt=1, so that the area under the pulse is given by A_(j). The observed signal consists of the detector output corrupted by noise, i.e.,

$\begin{matrix} {{s(t)} = {{r(t)} + {{w(t)}.}}} & (2) \end{matrix}$

The mathematical goal of pulse pile-up correction is to estimate ƒ_(A)(x) given a uniformlysampled, finite-length version of s(t). We assume throughout that the noise distribution w(t) is zero-mean Gaussian with known variance σ². We also assume the photon arrival times form a homogeneous Poisson process with a known rate. Consisting of the detector response R corrupted by a noise process W, where S_(k), R_(k) and W_(k) represent the kth elements of each series and where 0≤k<K. Let S, R and W be the uniformly sampled time-series corresponding to these signals,

$\begin{matrix} {S = {R + W}} & (3) \\ {= \left\{ {{s\left( t_{k} \right)}:{0 \leq k < K}} \right\}} & (4) \\ {R = \left\{ {{r\left( t_{k} \right)}:{0 \leq k < K}} \right\}} & (5) \\ {W = \left\{ {W_{t_{k}}:{0 \leq k < K}} \right\}} & (6) \\ {{{where}\mspace{14mu} t_{0}} < t_{1} < {.\;.\;.} < {t_{K - 1}.}} & \; \end{matrix}$

Summary of Pulse Processing Methodologies: Numerous approaches have been proposed over the decades to address the issue of pulse pile-up. Approaches can be broadly categorized into two types: time-domain based and energy-domain based. A popular strategy is to attempt to detect when pile-up has occurred in the time domain, and either reject or compensate for the affected pulses. Early spectroscopic systems adopted a rejectionbased approach along with matched filtering. The disadvantage of this approach is that an increasing proportion of pulses are rejected as the probability of pile-up grows. The system rapidly succumbs to paralysis, placing an upper limit on the count rate [1]. Strategies that compensate or correct for pile-up have grown in number with the increase of cheap computational power. These include template fitting [2], baseline subtraction [3], adaptive filtering [4, 5], sparse regression [6, 7] and more. These approaches all attempt to identify and compensate for pile-up in the time domain, and are generally best suited to systems with low pulse shape variation. The complexity of these approaches increases significantly with increasing variability between pulse shapes Φ_(j). It can be shown that any method that attempts to characterise individual pulses will suffer from pile-up. The best that these approaches can do is to reduce the onset of pile-up. Energy-based approaches attempt to address pile-up based on the statistics of an ensemble of pulses rather than individual pulses. They typically operate on histograms of estimated energy (the areas under clusters of pulses). The early work of Wielopolski and Gardner [8] and more recent extensions of their idea [9] operate primarily in the energy domain using ensemble-based strategies. Trigano et al. [10, 11] estimate the incident spectrum utilizing marginal densities from the joint distribution of the statistical properties of variable length clusters of pulses where the beginning and end of each cluster is detected. This circumvents the need to identify individual pulses, and is robust to pulse shape variation. Ilhe et al. [12] examine exponential shot-noise processes, restricting pulse shapes to a simple exponential to obtain tractable results. Further work [13] has been done to allow a wider range of pulse shapes. In both cases, knowledge is required of the pulse shape, along with estimates of the characteristic function and its derivative.

2 SUMMARY OF THE INVENTION

We chose an energy-based pile-up correction approach in order to i) avoid the limitations associated with the detection of individual pulses [14] and to ii) be able to handle pulse shape variation without undue increase in computational complexity. Rather than utilizing the joint distribution [10, 11] or shot-process [12] approaches, we recast the pile-up problem as a ‘decompounding’ problem of a compound Poisson process. A compound Poisson process is a discrete-time random process where each component consists of the sum of a random number of independent identically distributed random variables, where the number of random variables in each sum is Poisson distributed [15]. ‘Decompounding’ of a compound Poisson process is the task of using the random sums to estimate the distribution from which the random variables have been drawn. Buchmann and Grübel [16] formulated the decompounding of compound Poisson processes in the context of insurance claims and queuing theory. Decompounding of uniformly sampled compound Poisson processes has received some attention in recent times [16, 17, 18, 12, 19]. The context of these derivations frequently assume (reasonably) that each event is detectable (i.e., there is no ambiguity regarding the number of events), or that the density estimators are conditioned on at least one event occurring in each observation [20]. These assumptions are of limited value when addressing the spectroscopic pile-up problem.

The investigation of non-parametric decompounding without conditioning on event detection has received relatively little attention in the literature. Gugushvili [18] proposes a non-parametric, kernel-based estimator for the decompounding problem in the presence of Gaussian noise. In embodiments of the invention, the inventors have conceived that once a method for selecting the kernel bandwidth is obtained, along with a method for transforming the observed detector output to fit the mathematical model, this estimator can be readily extended and applied to a reformulation of the spectroscopic pile-up problem.

In accordance with a first broad aspect of the invention there is provided a method of determining a spectrum of energies of individual quanta of radiation received in a radiation detector, the method comprising the steps of: (1) obtaining a time series of digital observations from the radiation detector comprising pulses corresponding to the detection of the individual quanta; (2) computing spectrum sensitive statistics from the detector signal, the spectrum sensitive statistics defining a mapping from a density of amplitudes of the pulses to a spectrum sensitive statistics; (3) determining the spectrum by estimating the density of amplitudes of the pulses by applying an inversion of the mapping to the spectrum sensitive statistics.

In embodiments, the spectrum sensitive statistics may be based on a sum of the digital observations over a plurality of time intervals and the mapping may be defined using an approximate compound Poisson process, which may be augmented by a modelled noise. The mapping may be expressed as a relation between characteristic functions of the amplitudes, the spectrum sensitive statistics and the modelled noise. The characteristic functions of the spectrum sensitive statistics may be computed with the use of a histogram of the sum of the digital observations to which is applied an inverse Fourier transform. Computation of a characteristic function of the amplitudes may comprise the use of a low pass filter.

In a first embodiment, the plurality of time intervals are nonoverlapping and have have a constant length L and each interval is selected to encompass zero or more approximately entire clusters of the pulses. This may be accomplished by requiring a maximum value of the detector signal at a beginning and end of each time interval. In this embodiment, the compound Poisson process may be defined as a sum of the amplitudes of the pulses in each time interval. The mapping may be expressed as defined in equations (40) and (41), which may be augmented by windowing functions.

In a second embodiment, a plurality of intervals comprise a first set of nonoverlapping time intervals constant length L selected without regard to entirety of clusters of the pulses, and a 2nd set of nonoverlapping time intervals of constant length L1 less than L also selected without regard to entirety of clusters of the pulses. L is at least as long as a duration of the pulses and preferably L1 is less than the duration of the pulses. In this embodiment, the compound Poisson process may be defined as in Section 6. The mapping may be expressed as defined in Section 6. The second embodiment may utilise processes and calculations for each set of time intervals as defined for the set of time intervals in the first embodiment.

In embodiments, a data-driven strategy is used that results in a near optimal choice for a kernel parameter, which minimises the integrated-square-of-errors (ISE) of the estimated probability density function of incident photon energies.

According to a second broad aspect of the invention there is provided a method of estimating count rate of individual quanta of radiation received in a radiation detector, the method comprising the steps of: (1) obtaining a time series of digital observations from the radiation detector comprising pulses corresponding to the detection of the individual quanta; (2) computing spectrum sensitive statistics from the detector signal, the spectrum sensitive statistics using the intervals of constant length L and constant length L1 as described above in relation to the 1st broad aspect; (3) determining an estimate of a characteristic function of the compound Poisson process by formula (99); (4) estimating the count rate from the estimate of the characteristic function. Step (4) above may be achieved by using an optimisation routine or some other means to fit a curve, estimating a DC offset of a logarithm of the estimate of the characteristic function, or fitting a curve to the logarithm of the estimate of the characteristic function.

The rest of this application is organized as follows. Sections 3, 4 and 5 relate to the 1st embodiment of the 1st aspect of the invention. Sections 3 provides preliminary background, defines notation, outlines the mathematical model and gives a derivation of the estimator of the 1st embodiment, including modifications. Section 4 shows the performance of the modified estimator of the 1st embodiment for both simulated and experimental data, and discusses the results. Section 5 provides a conclusion for the 1st embodiment. Section 6 describe the 2nd embodiment, with reference to the 1st embodiment where relevant. Section 7 describes the 2nd aspect of the invention, a novel method of estimating count rate.

3 DERIVATION OF ESTIMATOR OF THE FIRST EMBODIMENT

The general approach we take to addressing pile-up is based on the following strategy; i) obtain statistics from s(t) that are sensitive to the distribution of incident photon energies, and estimate those statistics using the observed, finite-length sampled version of s(t), ii) obtain a mapping from the density of incident photon energies to the statistical properties of the observed statistics, iii) estimate the density of the incident photon energies by inverting the mapping. Section 3.1 describes our choice of statistics. Section 3.2 argues that these statistics (approximately) have the same distribution as a compound Poisson process. Section 3.3 introduces a decompounding technique for recovering the spectrum from these statistics. It is based on the decompounding algorithm in [18] but further developed to obtain near optimal performance in terms of the integrated square of error. Our approach to the pile-up problem follows the general theme of finding some statistics of s(t) that are sensitive to the underlying spectrum, estimating these statistics from a finite-time sampled version of s(t), then inverting the map that describes the statistical properties of these statistics given the underlying spectrum, thereby producing an estimate of the spectrum.

3.1 Choice of Statistic

We wish to obtain estimates of the photon energies from the observed signal given in (2). In typical modern spectroscopic systems, the detector output s(t) is uniformly sampled by an ADC. Without loss of generality, we assume the raw observations available to the algorithm are {s(k): k∈

_(≥0)}. Since identification of individual pulses can be difficult, we look instead for intervals of fixed length L∈

_(>0) containing zero or more clusters of pulses. Precisely, we define these intervals to be [T_(j), T_(j)+L) where

$\begin{matrix} {T_{0} = {\inf\left\{ {{k:\left| {s(k)} \middle| {\leq \epsilon} \right.},\left. {\left. {s\left( {k + L} \right)} \middle| {\leq \epsilon} \right.,{k \in {\mathbb{Z}}_{\geq 0}}} \right\}} \right.}} & (7) \\ {T_{j} = {\inf\left\{ {{{k:}❘\left. {s(k)} \middle| {\leq \epsilon} \right.},{\left. {\left. {s\left( {k + L} \right)} \middle| {\leq \epsilon} \right.,{{T_{j - 1} + L} \leq k},{k \in {\mathbb{Z}}_{\geq 0}}} \right\}.}} \right.}} & (8) \end{matrix}$

Here, ϵ is chosen as a trade-off between errors in the energy estimate and the probability of creating an interval. The value of ϵ should be sufficiently small to ensure the error in the estimate of total photon energy arriving within each interval is acceptably low, yet sufficiently large with respect to the noise variance to ensure a large number of intervals are obtained. Although the probability of partitioning the observed data into intervals approaches zero as the count-rate goes to infinity, this approach succumbs to paralysis at higher count-rates compared to pile-up rejection strategies based on individual pulses, since multiple photons are permitted to pile-up within in each interval. Section 4.2 describes the selection of L and ϵ for real data. Each interval contains an unknown, random number of pulses and may contain zero pulses.

We estimate the total photon energy x_(j) in the interval [T_(j), T_(j)+L) using the sampled raw observations. Since the area under each pulse is proportional to the photon energy A_(j) defined in (1), we let

$\begin{matrix} {x_{j} = {\sum_{k = T_{j}}^{T_{j} + L - 1}{s(k)}}} & (9) \end{matrix}$

The number of photon arrivals, the energy of each arriving photon and the detector output noise in each interval [T_(j), T_(j)+L) are assumed to be random and independent of other intervals. For pulse shapes with exponential decay, a small amount of the photon energy arriving in an interval may be recorded in the next interval. The amount of leakage is proportional to ϵ, and is negligible for sufficiently small ϵ. Consequently, the estimates x₁, x₂, . . . may be treated as the realization of a weakly-dependent, stationary process where each estimate is identically distributed according to the random variable X. This relationship is illustrated in FIG. 1 for the noise free case using a typical pulse shape.

3.2 Approximation with Compound Poisson Process

In this subsection we describe the distribution of X in terms of ƒ_(A)(x). We will then invert this in section 3.3, to obtain an estimator for the density ƒ_(A)(x). Using (9), (2), (1) and the fact that

(t) is causa we have

$\begin{matrix} {x_{j} = {{\sum_{\{{\ell:{\tau_{\ell} < T_{j}}}\}}{\sum_{k = T_{j}}^{T_{j} + L - 1}{a_{\ell}{\Phi_{\ell}\left( {k - \tau_{\ell}} \right)}}}} + {\sum_{\{{\ell:{T_{j} \leq \tau_{\ell} < {T_{j} + L}}}\}}{\sum_{k = T_{j}}^{T_{j} + L - 1}{a_{\ell}{\Phi_{\ell}\left( {k - \tau_{\ell}} \right)}}}} + {\sum_{k = T_{j}}^{T_{j} + L - 1}{w(k)}}}} & (10) \end{matrix}$

As justified below, this simplifies to

$\begin{matrix} {x_{j} \approx {y_{j} + z_{j}}} & (11) \\ {{{where}\mspace{11mu} y_{j}} = {\sum_{\{{\ell:{T_{j} \leq \tau_{\ell} < {T_{j} + L}}}\}}a_{\ell}}} & (12) \\ {{{and}\mspace{14mu} z_{j}} = {\sum_{k = T_{j}}^{T_{j} + L - 1}{{w(k)}.}}} & (13) \end{matrix}$

Both y_(j) and z_(j) are i.i.d. sequences of random variables. We denote their distributions by Y and Z. The distribution of Z is fully determined from the distribution of w(t), which is assumed zero-mean Gaussian with known variance σ². Moreover, Y is a compound Poisson process since the number of terms in the summation (number of photon arrivals in an interval of length L) has Poisson statistics. Equations (11)-(13) are justified as follows. The first term of (10) represents leakage from earlier intervals and is approximately zero. This is easily shown for Gaussian noise by performing a Taylor expansion about ϵ=0

$\begin{matrix} {{P{r\left( \left| {s(k)} \middle| {< \epsilon} \right. \right)}} = {{\frac{1}{2}{{erf}\left( \frac{{r(k)} + \epsilon}{\sqrt{2}\sigma} \right)}} - {\frac{1}{2}{{erf}\left( \frac{{r(k)} - \epsilon}{\sqrt{2}\sigma} \right)}}}} & (14) \\ {\approx {{0.79785\frac{\epsilon}{\sigma}e^{- \frac{{r{(k)}}^{2}}{2\sigma^{2}}}} + {{O\left( \epsilon^{3} \right)}.}}} & (15) \end{matrix}$

Thus there is a finite but small probability that some energy belonging to a previous interval will be included in the current estimate. In practice, this contribution is comparable to the noise for sufficiently small E. The third term in (10) is zero since

(t) is causal. The second term in (10) can be written as

$\begin{matrix} {{\sum\limits_{\{{\ell:{T_{j} \leq \tau_{\ell} < {T_{j} + L}}}\}}{a_{\ell}{\overset{T_{j} + L - 1}{\sum\limits_{k = T_{j}}}{\Phi_{\ell}\left( {k - \tau_{\ell}} \right)}}}} \approx {\sum\limits_{\{{\ell:{T_{j} \leq \tau_{\ell} < {T_{j} + L}}}\}}a_{\ell}}} & (16) \end{matrix}$

where we assume the pulse shapes

(t) are sufficiently smooth such that

Σ_(k = T_(j))^(T_(j) + L − 1)Φ_(ℓ)(k − τ_(ℓ)) ≈ ∫_(−∞)^(∞)Φ_(ℓ)(t)dt = 1.

It approximates the total energy of all the photons arriving in the interval [T_(j), T_(j)+L). Let ν_(j) designate the number of photon arrivals in the interval [T_(j), T_(j)+L). We assume ν_(j) is a realization of a homogeneous Poisson process with rate parameter λ, where λ is expressed in terms of the expected number of photons per interval of length L. Henceforth we shall assume that (11) holds exactly, and write

$\begin{matrix} {X = {Y + Z}} & (17) \end{matrix}$

Finally, we write x_(j) as

$\begin{matrix} {x_{j} \approx {y_{j} + z_{j}}} & (18) \end{matrix}$

where we assume Z has known variance σ². In this subsection we model the statistic of section 3.1 using a compound Poisson process. This allows us to derive an estimator for the density ƒ_(A)(x) in terms of observable quantities. The number of photons arriving in the interval [T_(j), T_(j)+L) is a Poisson random variable which we designate ν_(j). The total energy in the interval Y can be modelled as a compound Poisson process i.e.,

$\begin{matrix} {y_{j} = \left\{ \begin{matrix} {{\sum_{k = 0}^{v_{j} - 1}a_{({\ell_{j,1} + k})}},} & {\nu_{j} > 0} \\ 0 & {\nu_{j} = 0} \end{matrix} \right.} & (19) \\ {\nu_{j} \sim {P{n(\lambda)}}} & (20) \end{matrix}$

where

_(j,1)=min{

: T_(j)≤τ

<T_(j)+L} is the index of the first-photon arrival time in the interval, the arrival times are assumed ordered, and

representing photon energy are independent realizations of the random variable A with density function ƒ_(A)(x). The {ν_(j)} form a homogeneous Poisson process with rate parameter λ. The Poisson rate λ is expressed in terms of the expected number of photons per interval of length L.

$\begin{matrix} {y_{j} \approx {\sum_{k = T_{j}}^{T_{j} + L}{{r(k)}.}}} & (21) \end{matrix}$

The relationship between realizations of Y and the sampled detector response is illustrated in FIG. 1. The observed x_(j) can be approximated in terms of the

by substituting (2) into (9),

$\begin{matrix} {x_{j} = {\sum\limits_{k \in {\lbrack{{\hat{T}}_{j},{{\hat{T}}_{j} + L}}\rbrack}}\left( {{r(k)} + {w(k)}} \right)}} & (22) \\ {x_{j} = {{\sum_{k = T_{j}}^{T_{j} + L}{r(k)}} + {\sum_{k = T_{j}}^{T_{j} + L}{w(k)}}}} & (23) \\ {\approx {y_{j} + z_{j}}} & (24) \end{matrix}$

where z_(j) is the realization of the unobservable random variable Y that represents the photon energy in an interval of the discrete-time detector response,

$\begin{matrix} {Y_{j} = {\sum\limits_{k \in {\lbrack{T_{j},{T_{j} + L}}\rbrack}}R_{k}}} & (25) \end{matrix}$

where z_(j) is a realization of Z, an independent random variable representing errors in the sampling process and estimation of T_(j). We assume Z has known variance σ². With these definitions of X and Y, the number of intervals which can be found in a finite length of detector output is a random variable N. At high count-rates this approach succumbs to paralysis, as the probability of being able to partition the observed data into intervals approaches zero. The onset of paralysis occurs at higher count-rates compared to pile-up rejection based strategies, since multiple photons are permitted to pile-up within in each interval. Assume the time-series defined in (3)-(6) has been sampled uniformly. Without loss of generality, assume unit sample intervals beginning at t₀=0 i.e., t_(k)=k, 0≤k<K. Let R be a discrete-time random process representing the sampled detector response of (1). Let Y={Y_(j): 0≤j<N} be a discrete-time random process whose components Y_(j) represent the total photon energy arriving during a fixed time interval. A compound Poisson process can be used to model Y, i.e.,

$\begin{matrix} {Y_{j} = \left\{ \begin{matrix} {{\sum_{k = 1}^{v_{j}}A_{k}},} & {\nu_{j} > 0} \\ 0 & {\nu_{j} = 0} \end{matrix} \right.} & (26) \\ {\nu_{j} \sim {P{n(\lambda)}}} & (27) \end{matrix}$

where ν_(j) is an independent Poisson random variable, and A_(k) are independent identically distributed random variables with density function ƒ_(A)(x). The {ν_(j)} form a homogeneous Poisson process with rate parameter λ. The process Y is not directly observable. Assume the pulse shape Φ(t) has finite support. Let

_(A)(t) be the indicator function for the set A. Let the pulse length

be given by

=sup({t: Φ>0})−inf({t: Φ(t)>0}). Let S=R+W={s(k):0≤k<K} be a discrete-time random process representing the observed detector output given by (2). It consists of the detector response R corrupted by a noise process W. Without loss of generality, we assume unit sample intervals. From the observations S we form the process X, where

$\begin{matrix} {X_{j}\overset{\bigtriangleup}{=}{Y_{j} + Z_{j}}} & (28) \end{matrix}$

and where Z_(j) is a random variable from an independent noise process of known variance σ². A simple model for testing theory is obtained when we let the pulse shape Φ(t)=

_((0,1))(t) in (1), in which case we let X_(j)=S_(j), and N is simply the sample length K. Obtaining X_(j) from S is more complicated for real data. In that case we partition the process S into non-overlapping blocks of length L, where L>

. The Poisson rate λ is expressed in photons per block. The start of each block T_(j)∈

is chosen such that the total energy of any pulse is fully contained within the block in which it arrives

$\begin{matrix} {T_{j} = {\min\mspace{11mu}\left\{ {{{k\ \text{:}\ R_{k}} = 0},\ {R_{k + L} = 0},\ {{T_{j - 1} + L} < k < {K - L}}} \right\}}} & (29) \end{matrix}$

FIG. 1 shows that

$Y_{j} = {\sum\limits_{k = T_{j}}^{T_{j} + L}\;{R_{k}.}}$

We let

$X_{j} = {\sum\limits_{k = {\hat{T}}_{j}}^{{\hat{T}}_{j} + L}S_{k}}$

$\begin{matrix} {X_{j} = {\sum\limits_{k = {\hat{T}}_{j}}^{{\hat{T}}_{j} + L}S_{k}}} & (30) \end{matrix}$

where {circumflex over (T)}_(j) is an estimate of T_(j). Section 4.2 describes the selection of L and ϵ for real data. With this definition of X_(j), the number of components in Y becomes a random variable for a given sample length K. At high count-rates this approach succumbs to paralysis, as the probability of being able to create a block approaches zero. The onset of paralysis occurs at higher count-rates compared to pile-up rejection based strategies, since multiple photons are permitted to pile-up within in each block. Let Y={Y_(j): 0≤j<N} be a discrete-time random process whose components Y_(j) are given by

$\begin{matrix} {Y_{j} = {\sum\limits_{k = T_{j}}^{T_{j} + L}R_{k}}} & (31) \\ {T_{j} = {\min\limits_{k > {T_{j - 1} + L}}\left\{ {{R_{k} < d},\ {R_{k + L} < d}} \right\}}} & (32) \end{matrix}$

where L∈

is a constant chosen such that h and d is a small threshold value close to zero. The random variable Y_(j) thus represents the total photon energy arriving during a fixed time interval of length L. The value of d ensures the signal associated with photon arrivals is very small at the start and end of each interval. This is illustrated in FIG. 1. A compound Poisson process can be used to model Y, i.e.,

$\begin{matrix} {Y_{j} = \left\{ \begin{matrix} {{\sum\limits_{k = 1}^{v_{j}}A_{k}},} & {v_{j} > 0} \\ 0 & {v_{j} = 0} \end{matrix} \right.} & (33) \\ {v_{\lambda} = \left\{ {v_{j}:{0 \leq j < N}} \right\}} & (34) \\ {v_{j} \sim {{Pn}(\lambda)}} & (35) \end{matrix}$

where ν_(λ) is a homogeneous Poisson process with rate parameter λ, and A_(k) are independent identically distributed random variables with density function ƒ_(A)(x). Let S=R+W be a discrete-time random process representing the sampled detector output given by (2). It consists of the detector response R corrupted by a noise process W. The process Y is not directly observable. Using (2), (25) and (32), we model observations by the process X={X_(j): 0≤j<N}, i.e.,

$\begin{matrix} {{X_{j} = {\sum\limits_{k = T_{j}}^{T_{j} + L}S_{k}}}\mspace{565mu}} & {(36)} \\ {= {\sum\limits_{k = T_{j}}^{T_{j} + L}\left( {R_{k} + W_{k}} \right)}} & {(37)} \\ {= {Y_{j} + {\sum\limits_{k = T_{j}}^{T_{j} + L}W_{k}}}} & {(38)} \\ {\overset{\Delta}{=}{Y_{j} + Z_{j}}} & {(39)} \end{matrix}$

-   -   where Z is a noise process of known variance σ². All the random         variables (ν_(j), A₁, . . . , A_(ν) _(j) , Z_(j)) involved in         modelling a given observation X_(j) are assumed independent. Let         X₁, X₂, . . . , X_(N) be N independent, identically distributed         observations. Let X, Y, Z, A be the collections of X_(j), Y_(j),         Z_(j), A_(j): 0≤j<N. Let the corresponding characteristic         functions be ϕ_(X), ϕ_(Y), ϕ_(Z), ϕ_(A).

3.3 Basic Form of Estimator

We seek to invert the mapping from the distribution of photon energy A to the distribution of X. Our strategy is to first obtain the characteristic function of X in terms of ƒ_(A), then invert the mapping assuming the count-rate and noise characteristics are known. Let ϕ_(X), ϕ_(Y), ϕ_(Z), ϕ_(A) be the characteristic functions of X, Y, Z, A. It is well known [15] that for the compound Poisson process Y with rate λ,

$\begin{matrix} {{\phi_{Y}(u)} = {e^{- \lambda}e^{\lambda_{\phi_{A}{(u)}}}}} & (40) \end{matrix}$

and since X=Y+Z then

$\begin{matrix} {{\phi_{X}(u)} = {{\phi_{Y}(u)}{{\phi_{Z}(u)}.}}} & (41) \end{matrix}$

Given the observations x_(j) we can form an empirical estimate {circumflex over (ϕ)}_(X) of the characteristic function of X. Treating this as the true characteristic function, we can invert (40), (41) to obtain the characteristic function of A and then take the Fourier transform to find the amplitude spectrum ƒ_(A). Specifically, using (40), (41) and exploiting the assumption that Z is Gaussian to ensure ϕ_(Z)(u) will be non-zero ∀u∈

, we let γ:

→

be the curve described by

$\begin{matrix} {{\gamma(u)} = \frac{\phi_{X}(u)}{e^{- \lambda}{\phi_{Z}(u)}}} & (42) \\ {= e^{{\lambda\phi}_{A}{(u)}}} & (43) \end{matrix}$

Temporarily assuming ∀u, γ(u)≠0, after taking the distinguished logarithm of (43) and rearranging we have

$\begin{matrix} {{\phi_{A}(u)} = {\frac{1}{\lambda}{{dlog}(\gamma)}\mspace{11mu}{(u).}}} & (44) \end{matrix}$

Ideally, ƒ_(A) is recovered by taking a Fourier transform

$\begin{matrix} {{f_{A}(x)} = {\int_{- \infty}^{\infty}{e^{{- i}\; 2\pi\;{ux}}{\phi_{A}(u)}\mspace{11mu}{du}}}} & (45) \end{matrix}$

The basic form of our proposed estimator is given in (88) and is derived from (45) via a sequence of steps. First, ϕ_(X) is estimated from the data (Step 1). Simply substituting this estimate for ϕ_(X) in (42) does not produce an ISE optimal estimate of γ. The approximate ISE is obtained from an approximate estimate of the error distribution of ϕ_(X) (Step 2). We then determine a sensible windowing function G(u) (in Step 3) and estimate γ by

$\begin{matrix} {{{\overset{\hat{}}{\gamma}(u)} = {{G(u)}\frac{{\overset{\hat{}}{\phi}}_{X}(u)}{e^{- \lambda}{\phi_{Z}(u)}}}}.} & (46) \end{matrix}$

The windowing function G(u) is designed to minimise the approximate ISE between ƒ_(A) and our estimate of ƒ_(A) based on (44), (45) and (46), but with γ in (44) replaced by (46). A similar idea is used for estimating ϕ_(A) from (44): a weighting function H(u) is found (in Step 4) such that replacing ϕ_(A) in (45) by

$\begin{matrix} {{{\hat{\phi}}_{A}(u)} = {{H(u)}\frac{1}{\lambda}{dlog}\mspace{11mu}\left( \overset{\hat{}}{\gamma} \right)}} & (47) \end{matrix}$

produces a better estimate of ƒ_(A) than using the unweighted estimate 1/λd log ({circumflex over (γ)}). Finally, the weighting function H(u) is modified (in Step 5) to account for the integral in (45) having to be replaced by a finite sum in practice. The following subsections expand on these five steps.

3.4 Estimating ϕ_(X)

An estimate of ϕ_(X)(u) is required to estimate γ(u). In this subsection we define a histogram model and describe our estimation of ϕ_(X)(u) based on a histogram of the x_(j) values. Assume N intervals (and corresponding x_(j) values) have been obtained from a finite length data sample. Although the empirical characteristic function

$\begin{matrix} {{{\overset{\hat{}}{\phi}}_{X\mspace{11mu}{emp}}(u)} = {\frac{1}{N}{\sum\limits_{j = 0}^{N - 1}e^{iux_{j}}}}} & (48) \end{matrix}$

provides a consistent, asymptotically normal estimator of the characteristic function [21], it has the disadvantage of rapid growth in computational burden as the number of data points N and the required number of evaluation points u∈

increases. Instead, we use a histogram based estimator that has a lower computational burden. Assume that a histogram of the observed X values is represented by the 2M×1 vector n, where the count in the mth bin is given by

n m = ∑ k = 0 N - 1 ⁢ [ m - 0.5 , m + 0.5 ) ⁢ ( x k ) , m ∈ { - M , … ⁢ , M - 1 } . ( 49 )

All bins of the histogram have equal width. The bin-width is chosen in relation to the magnitude of the x_(j) values. Since the effect of choosing a different bin width is simply equivalent to scaling the x_(j) values, we assume the bin-width to be unity without loss of generality. The bins are apportioned equally between non-negative and negative data values. The number of histogram bins 2M influences the estimator in various ways, as discussed in later subsections. For now, it is sufficient to assume that 2M is large enough to ensure the histogram includes all x_(j) values. We estimate ϕ_(X)(u) by forming a histogram of scaled x_(j) values and take the inverse Discrete Fourier transform i.e.,

$\begin{matrix} {{\phi_{X}(u)} = {\sum\limits_{m = {- M}}^{M - 1}{\frac{n_{m}}{N}{e^{\frac{i\; 2\pi\mspace{11mu} u\mspace{11mu} m}{2M}}.}}}} & (50) \end{matrix}$

This is a close approximation of the empirical characteristic function but where x_(j) terms have been rounded to the nearest histogram bin centre (and u contracted by a factor of 2π). The term n_(m) simply counts the number of rounded terms with the same value. Clearly, this function can be efficiently evaluated at certain discrete points u∈−M, . . . , M−1 using the fast Fourier Transform (FFT).

3.5 Error Distribution of {circumflex over (ϕ)}_(X)

The design of the filters G(u) and H(u) in (46) and (47) rely on the statistics of the errors between {circumflex over (ϕ)}_(X) and the true characteristic function. In this subsection we define and describe the characteristics of these errors. We assume the density function ƒ_(X) is sufficiently smooth (i.e., |d^(n)ƒ_(X)(u)/du^(n)|≤C_(n)∈

∀n∈

) and that the width of the histogram bins are sufficiently small (relative to the standard deviation of the additive noise Z) such that the errors introduced by rounding x_(j) values to the centre of each histogram bin are approximately uniformly distributed across each bin, have zero mean and are small relative to the peak spreading caused by Z. In other words, the source of error arising from the binning of x_(j) values is considered negligible. Due to both the statistical nature of Poisson counting and the expected count in each bin being non-integer (

[n_(m)]∈

_(≥0)), discrepancies exist between the observed number of counts in any given histogram bin and the expected number of counts for that bin. We combine these two sources of error in our model and refer to it as ‘histogram noise’. We emphasize that this noise is distinct from the additive noise Z modelled in (11), which causes peak spreading in the histogram. Let the probability that a realization of X falls in the m-th bin be

$\begin{matrix} {{px}_{m} = {\Pr\left( {{m - 0.5} \leq X < {m + 0.5}} \right)}} & (51) \end{matrix}$

Let the normalized histogram error ϵ_(m) in the m-th bin be the difference between the observed count n_(m) and the expected count

[n_(m)]=Npx_(m) in the mth bin, relative to the total counts in the histogram N i.e.,

$\begin{matrix} {\epsilon_{m} = \frac{n_{m} - {Npx_{m}}}{N}} & (52) \end{matrix}$

Using (50), (51) and (52) we have

$\begin{matrix} \begin{matrix} {{{\overset{\hat{}}{\phi}}_{X}(u)} = {\sum_{m = {- M}}^{M - 1}{\frac{n_{m}}{N}e^{i\frac{2\pi um}{2M}}}}} \\ {= {{\sum_{m = {- M}}^{M - 1}{{px}_{m}e^{i\frac{2\pi um}{2M}}}} + {\sum_{m = {- M}}^{M - 1}{\epsilon_{m}e^{i\frac{2\pi um}{2M}}}}}} \\ {\approx {{\phi_{X}(u)} + {\phi_{\epsilon}(u)}}} \end{matrix} & \begin{matrix} \begin{matrix} \begin{matrix} (53) \\ (54) \end{matrix} \\ \; \end{matrix} \\ (55) \end{matrix} \end{matrix}$

If the histogram is modelled as a Poisson vector, show that

$\begin{matrix} {{{\mathbb{E}}\left\lbrack \epsilon_{i} \right\rbrack} = 0} & (56) \\ {{{\mathbb{E}}\left\lbrack {\epsilon_{i}\epsilon_{j}} \right\rbrack} = \left\{ \begin{matrix} \frac{{px}_{j}}{N} & {i = j} \\ 0 & {i \neq j} \end{matrix} \right.} & (57) \\ {{\mathbb{E}}\left\lbrack {\left. \left. \phi_{\epsilon} \right|^{2} \right\rbrack = {\frac{1}{N}.}} \right.} & (58) \end{matrix}$

Since the characteristics of the histogram noise can be expressed in terms of the total number of observed intervals N, the impact of using observation data of finite length may be accounted for by incorporating this information into the design of G(u) and H(u).

3.6 Estimating γ

Having obtained {circumflex over (ϕ)}_(X), the next task is to estimate γ. Rather than substitute {circumflex over (ϕ)}_(X)(u) for ϕ_(X)(u) in (42), we instead use (46) as the estimator, which requires us to choose a windowing function G(u). In this subsection we attempt to find a function G(u) that is close to optimal. When the distribution of errors in {circumflex over (ϕ)}_(X)(u) are considered, the windowing function G(u)=G_(opt)(u) that results in the lowest ISE estimator of the form given in (46) is

$\begin{matrix} {{G_{opt}(u)} = \frac{1}{\sqrt{1 + \frac{e^{{- 2}{\lambda\Re}{\{{\phi_{A}{(u)}}\}}}}{Ne^{{- 2}\lambda}{{\phi_{Z}(u)}}^{2}}}}} & (59) \end{matrix}$

where

{z} denotes the real component of z∈

. We cannot calculate G_(opt)(u) since ϕ_(A)(u) is unknown, so instead we attempt to find an approximation. We let

$\begin{matrix} {{G(u)} = {\frac{1}{\sqrt{1 + \frac{1}{Ne^{{- 2}\lambda}{\phi_{Z}}^{2}(u)}}}.}} & (60) \end{matrix}$

This is justified by considering the magnitude of the relative error between the functions g_(opt)(u) and g₁(u) where

$\begin{matrix} {{g_{opt}(u)} = {1 + \frac{e^{{- 2}{\lambda\Re}{\{{\phi_{A}{(u)}}\}}}}{Ne^{{- 2}\lambda}{\phi_{Z}}^{2}(u)}}} & (61) \\ {{g_{1}(u)} = {1 + {\frac{1}{\left. {Ne^{{- 2}\lambda}} \middle| \phi_{Z} \middle| {}_{2}(u) \right.}.}}} & (62) \end{matrix}$

The magnitude of the relative error is given by

$\begin{matrix} {{\frac{g_{opt} - g_{1}}{g_{1}}} = {{\frac{e^{{- 2}{\lambda\Re}{\{\phi_{A}\}}} - 1}{\left. {Ne^{{- 2}\lambda}} \middle| \phi_{Z} \middle| {}_{2}{+ 1} \right.}}.}} & (63) \end{matrix}$

Since

{ϕ_(A)}∈[−1,1], we see the right hand side of (63) is maximum when

{ϕ_(A)(u)}=−1. The relative error is thus bound by

$\begin{matrix} {{\frac{g_{opt} - g_{1}}{g_{1}}} \leq \frac{e^{2\lambda} - 1}{\left. {Ne^{{- 2}\lambda}} \middle| \phi_{Z} \middle| {}_{2}{+ 1} \right.}} & (64) \end{matrix}$

which justifies the approximation when λ is small, or when N|ϕ_(Z)|²(u)>>e^(4λ). Furthermore, we note that the above bound is quite conservative. The distribution of photon energies in spectroscopic systems can typically be modelled as a sum of K Gaussian peaks, where the kth peak has location μ_(k) and scale σ_(k) i.e.,

$\begin{matrix} {{{f_{A}(x)} = {\overset{K - 1}{\sum\limits_{k = 0}}{\alpha_{k}\frac{1}{\sqrt{2\pi}\sigma_{k}}e^{- \frac{{({x - \mu_{k}})}^{2}}{2\sigma_{k}^{2}}}}}}{where}} & (65) \\ {{\sum_{k = 0}^{K - 1}\alpha_{k}} = 1.} & (66) \end{matrix}$

Consequently, the characteristic function will have the form

$\begin{matrix} {{\phi_{A}(u)} = {\sum_{k = 0}^{K - 1}{\alpha_{k}e^{{- 2}\pi^{2}\sigma_{k}^{2}u^{2}}{e^{i2\pi\mu_{k}u}.}}}} & (67) \end{matrix}$

i.e., oscillations within an envelope that decays as e^(−cu) ² for somec>0. The upper bound given by (64) is quite conservative since |

{ϕ_(A)}|>>1 for most values of u. The approximation error will be significantly smaller at most evaluation points across the spectrum. Having chosen G(u), we can form an estimate of γ using (46). The windowing function reduces the impact of histogram noise arising from the finite number of data samples. For large values of Ne^(−2λ)|ϕ_(Z)(u)|², the impact of windowing is negligible and the estimator is essentially the same as using (42) directly. However, in the regions where

$\begin{matrix} {{\ln N} < {{2\lambda} + {4\pi^{2}\sigma^{2}u^{2}}}} & (68) \end{matrix}$

the windowing becomes significant, and acts to bound our estimate of γ i.e., Using the fact that the noise Z is Gaussian (so ϕ_(Z)(u)∈

and hence |ϕ_(Z)|²=ϕ_(Z) ²), and since e^(−2λ)>0 we see that

$\begin{matrix} {\begin{matrix} {{{\overset{\hat{}}{\gamma}(u)}} = \left| {\frac{{\overset{\hat{}}{\phi}}_{X}(u)}{e^{- \lambda}{\phi_{Z}(u)}}\frac{1}{\sqrt{1 + \frac{1}{Ne^{{- 2}\lambda}{\phi_{Z}}^{2}(u)}}}} \right|} \\ {= {{{{{\overset{\hat{}}{\phi}}_{X}(u)}\frac{1}{\sqrt{{e^{{- 2}\lambda}{\phi_{Z}^{2}(u)}} + \left( \frac{\phi_{Z}^{2}(u)}{N{\phi_{Z}}^{2}(u)} \right)}}}} <}} \\ {\sqrt{N}} \end{matrix}.} & \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} (69) \\ \; \end{matrix} \\ \; \end{matrix} \\ \; \end{matrix} \\ (70) \end{matrix} \\ \; \end{matrix} \\ \; \end{matrix} \\ (71) \end{matrix} \end{matrix}$

This ensures the argument to the distinguished logarithm in (47) remains finite even though lim_(u→∞)ϕ_(Z)(u)=0.

3.7 Estimating ϕ_(A)

Once {circumflex over (γ)} has been obtained, we proceed to estimate ϕ_(A) using (47). This requires another windowing function H(u). In this subsection we find a function H(u) for estimating ϕ_(A) that is close to ISE optimal. We begin by defining a function ψ(u) for notational convenience

$\begin{matrix} {{\psi(u)} = {\frac{1}{G(u)}e^{- \lambda}{{\phi_{Z}(u)}.}}} & (72) \end{matrix}$

The ISE is minimized when H(u)=H_(opt)(u), where the optimal filter H_(opt)(u) is given by

$\begin{matrix} \begin{matrix} {{H_{opt}(u)} = {\frac{\phi_{A}(u)}{{\overset{\hat{}}{\phi}}_{A}(u)} = \frac{\frac{1}{\lambda}d{\log\left( \frac{\phi_{X}}{e^{- \lambda}\phi_{Z}} \right)}(u)}{\frac{1}{\lambda}d{\log\left( \frac{\phi_{X} + \phi_{\epsilon}}{\psi} \right)}(u)}}} \\ {= \frac{\phi_{A}(u)}{{\phi_{A}(u)} + {\frac{1}{\lambda}d{\log\left( \frac{\phi_{X}}{\phi_{X}} \right)}(u)} + {\frac{1}{\lambda}d{\log(G)}(u)}}} \end{matrix} & \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} (73) \\ \; \end{matrix} \\ \; \end{matrix} \\ \; \end{matrix} \\ (74) \end{matrix} \\ \; \end{matrix} \end{matrix}$

Again, we cannot calculate the optimal filter by using (73)-(74) since ϕ_(X)(u), ϕ_(A)(u) and ϕ_(ϵ)(u) are unknown. We instead make the following observations to obtain an approximation of the ISE-optimal filter.

3.7.1 Initial Observations

The optimal filter remains close to unity as long as the estimated {circumflex over (ϕ)}_(A)(u) remains close to the true value of ϕ_(A)(u). This will invariably be the case for small values of u since

$\begin{matrix} \begin{matrix} {{{\mathbb{E}}\left\lbrack {{\phi_{\epsilon}(u)}} \right\rbrack} \approx \sqrt{\frac{\pi}{4N}} ⪡} \\ {{{\phi_{X}(u)}} \approx {1\mspace{14mu}{for}\mspace{14mu}{small}\mspace{14mu} u}} \end{matrix} & \begin{matrix} \begin{matrix} (75) \\ \; \end{matrix} \\ (76) \end{matrix} \end{matrix}$

Furthermore, equation (73) shows that if |ϕ_(ϵ)(u)|≤≤|ϕ_(X)(u)|, then {circumflex over (ϕ)}_(X)(u)=ϕ_(X)(u)+ϕ_(ϵ)(u)≈ϕ_(X)(u) so H_(opt)(u)≈1. For larger values of u, when the magnitude of |ϕ_(X)(u)| becomes comparable to or less than |ϕ_(ϵ)(u)|, the estimator

${{\overset{\hat{}}{\phi}}_{A}(u)} = {\frac{1}{\lambda}d\;{\log\left( {\left( {\phi_{X} + \phi_{\epsilon}} \right)/\psi} \right)}(u)}$

is dominated by noise and no longer provides useful estimates of ϕ_(A)(u). In the extreme case |ϕ_(X)(u)|<<ϕ_(ϵ)(u)| so |{circumflex over (ϕ)}_(X)(u)|≈|ϕ_(ϵ)(u)| and hence

$\begin{matrix} {{\mathcal{R}\left\{ {\hat{\phi}}_{A} \right\}} \approx {\frac{1}{\lambda}\ln{\frac{\phi_{\epsilon}}{\psi}}}} & (77) \end{matrix}$

The window H(u) should exclude these regions from the estimate, as the bias introduced in doing so will be less than the variance of the unfiltered noise. Unfortunately the estimate of ϕ_(A)(u) can be severely degraded well before this boundary condition is reached, so (77) is not particularly helpful. A more useful method for detecting when noise begins to dominate is as follows.

3.7.2 Filter Design Function

Further manipulation of (67) shows that for typical spectroscopic systems, the magnitude of ϕ_(A) will have the form

$\begin{matrix} {{{\phi_{A}}^{2}(u)} = {{\sum\limits_{k = 0}^{K - 1}\;{\alpha_{k}^{2}e^{{- 4}\pi^{2}\sigma_{k}^{2}u^{2}}}} + {\underset{j \neq k}{\sum\limits_{k = 0}^{K - 1}\;\sum\limits_{j = 0}^{K - 1}}\;\alpha_{k}\alpha_{j}\mspace{14mu}\cos\mspace{14mu}\left( {2{\pi\left( {\mu_{k} - \mu_{j}} \right)}u} \right)e^{{- 2}{\pi^{2}{({\sigma_{k}^{2} + \alpha_{j}^{2}})}}u^{2}}}}} & (78) \end{matrix}$

i.e.; a mean component that decays according to the peak widths σ_(k), and a more rapidly decaying oscillatory component that varies according to the location of the spectral peaks μ_(k). In designing the window H(u), we are interested in attenuating the regions in |{circumflex over (ϕ)}_(A)| where |ϕ_(A)|²≲|ϕ_(ϵ)/ψ|², i.e., where the signal power is less than the histogram noise that has been enhanced by the removal of ϕ_(Z) during the estimation of γ. To obtain an estimate of |ϕ_(A)|, a low-pass, Gaussian shaped filter H_(lpf)(u) is convolved with |{circumflex over (ϕ)}_(A)| to attenuate all but the slowly varying, large scale features of |{circumflex over (ϕ)}_(A)|. We denote this |{circumflex over (ϕ)}_(Asmooth)|(u)

$\begin{matrix} {{{{\hat{\phi}}_{Asmooth}}(u)} = {{{\frac{1}{\lambda}d\;{\log\left( \hat{\gamma} \right)}(u)}}\mspace{14mu}\bigstar\mspace{14mu}{{H_{1{pf}}(u)}.}}} & (79) \end{matrix}$

We see that |ϕ_(ϵ)(u)| has a Rayleigh distribution with scale parameter

$\sigma_{Ray} = {\frac{1}{\sqrt{2N}}.}$

Consequently

$\begin{matrix} {{\frac{1}{\lambda}\frac{{\phi_{\epsilon}(u)}}{\psi(u)}} \sim {{{Rayleigh}\left( {\sigma_{Ray} = \frac{1}{{{\lambda\psi}(u)}\sqrt{2N}}} \right)}.}} & (80) \end{matrix}$

It is well known that the cumulative distribution function of a Rayleigh distributed random variable X_(Ray) is given by

$\begin{matrix} {{F_{Ray}\left( {x;\sigma_{Ray}} \right)} = {\Pr\left( {{X_{Ray} < x};\sigma_{Ray}} \right)}} & (81) \\ {= {1 - {e^{- \frac{x^{2}}{2\sigma_{Ray}^{2}}}.}}} & (82) \end{matrix}$

Hence, to assist with computing the window H(u) we will make use of the function

$\begin{matrix} {{\alpha_{\min}(u)} = {1 - e^{{- \frac{1}{2}}{{\hat{\phi}}_{Asmooth}}{(u)}^{2}\lambda^{2}{\psi{(u)}}^{2}2N}}} & (83) \\ {\approx {\Pr\left( {{\frac{1}{\lambda}\frac{{\phi_{\epsilon}(u)}}{\psi(u)}} < {{{\hat{\phi}}_{Asmooth}}(u)}} \right)}} & (84) \end{matrix}$

to control the shape of H(u). The function α_(min)(u) provides an indication of how confident we can be that the estimate {circumflex over (ϕ)}_(A)(u) contains more signal energy than noise energy. The approximation in (84) arises from the fact that |{circumflex over (ϕ)}_(Asmooth)| is also a random variable slightly affected by the noise ϵ. On occasion—particularly for larger values of |u|—the histogram noise may result in sufficiently large values of α_(min)(u) to give a false sense of confidence, and potentially allow noisy results to corrupt the estimate of ϕ_(A). To overcome this problem, the function was modified to be uni-modal in u

$\begin{matrix} {{\alpha_{mod}(u)} = {\inf\left\{ {{\alpha_{\min}(\upsilon)},\ {{\upsilon } < {u}}} \right\}}} & (85) \end{matrix}$

This modification was justified on the assumption that Gaussian noise causes ϕ_(Z)(u) to be decreasing in |u|. Consequently we expect

[|ϕ_(ϵ)(u)|/ψ(u)] to be increasing in |u|. If we ignore the local oscillations in ϕ_(A)(u) that are due to peak locations in ƒ_(A)(x), the envelope approximated by the smoothed |ϕ_(Asmooth)|(u) will be non-increasing in |u|. Equation (74) indicates the optimal window has the form λϕ_(A)(u)/(λϕ_(A)(u)+d log({circumflex over (ϕ)}_(X)/ϕ_(X))(u)+d log(G)(u), so the overall window shape will be decreasing in |u|. Hence, if the estimated characteristic function in the region of some u₀ (where the signal to noise ratio is high) has determined that the window value should be H(u₀)<1, then it is reasonable to reject the suggestion that in the region u₁>u₀ (where the signal to noise ratio will be worse) that H(u₁)>H(u₀). Using the knowledge that |H_(opt)(u)| should be close to unity for small |u|, close to zero for large |u|, and should ‘roll off’ as the signal-to-noise-ratio decreases—we consider two potential windowing functions as approximations of H_(opt)(u).

3.7.3 Rectangle Window

The indicator function provides a very simple windowing function

$\begin{matrix} {{H(u)} = {1_{\{{{\alpha_{mod}{(u)}} > \alpha_{0}}\}}{(u).}}} & (86) \end{matrix}$

The threshold value α₀ determines the point at which cut-off occurs, and can be selected manually as desired (e.g., α₀=0.95). Once the threshold is chosen, the estimator exhibits similar ISE performance regardless of peak locations in the incident spectra. Rather than requiring the user to select a window width depending on the incident spectrum¹, the width of the window is automatically selected by the data via α_(mod)(u). While simplicity is the primary advantage of the rectangular window, the abrupt transition region provides a poor model for the roll-off region of the optimal filter. The second filter shape attempts to improve on that. ¹Gugushvili [18] proposed a non-parametric estimator for the general decompounding problem in which a rectangular windowing scheme was used. It requires manual selection of window width, which varies with ϕ_(A).

3.7.4 Logistic Window

A window based on the logistic function attempts to model smoother roll-off. It is given by

$\begin{matrix} {{H(u)} = \frac{1 + e^{- {\beta_{0}{({1.0 - \alpha_{0}})}}}}{1 + e^{- {\beta_{0}{({{\alpha_{mod}{(u)}} - \alpha_{0}})}}}}} & (87) \end{matrix}$

where α₀ again acts as a threshold of acceptance of the hypothesis that the signal energy is greater than the noise energy in the estimate {circumflex over (ϕ)}_(A)(u). The rate of filter roll off in the vicinity of the threshold region is controlled by β₀>0. This provides a smoother transition region than the rectangle window, reducing Gibbs oscillations in the final estimate of ϕ_(A). Once again, although the parameters α₀, β are chosen manually, they are much less dependent on ϕ_(A) and can be used to provide close to optimal filtering for a wide variety of incident spectra. Typical values used were α₀=0.95, β₀=40.0. The performance of the rectangle and logistic window functions are compared in section 4.

3.8 Estimating ƒ_(A)

Having designed a window function H(u) and thus an estimator {circumflex over (ϕ)}_(A)(u), the final task is to estimate ƒ_(A)(x) by inverting the Fourier transform. This sub-section describes several issues that arise with numerical implementation. Firstly, it is infeasible to evaluate {circumflex over (ϕ)}_(X), {circumflex over (γ)}(u) and {circumflex over (ϕ)}_(A) numerically on the whole real line. Instead we estimate it at discrete points over a finite interval. The finite interval is chosen sufficiently large such that a tolerably small error is incurred as a result of excluding signal values outside the interval. This is justified for ƒ_(A)(x) being a Gaussian mixture, since the magnitudes of ϕ_(X) and ϕ_(A) will decay as e^(−cu) ² for some c>0. The Fast Fourier Transform (FFT) is used to evaluate {circumflex over (ϕ)}_(X) at discrete points, and hence also determines the points where {circumflex over (γ)}(u) and {circumflex over (ϕ)}_(A) are evaluated. Likewise, the FFT is used to evaluate the final estimate of {circumflex over (ƒ)}_(A) at discrete points. In order to use the FFT, the signals outside the interval should be sufficiently small to reduce impact of aliasing. The evaluation points also need to be sufficiently dense to avoid any ‘phase wrap’ ambiguity when evaluating d log({circumflex over (γ)})(u). Both these objectives can be achieved by increasing the number of bins 2M in the histogram (zero-padding) until a sufficiently large number of bins is attained. As M increases, the sampling density of {circumflex over (γ)} increases, which allows phase wrapping to be detected and managed. A larger M also allows aliasing (caused by the Gaussian shaped tails of |ϕ_(X)|) to be negligible. Typically a value of M was chosen as the smallest power of two sufficiently large such that the non-zero values of the histogram were confined to the ‘lower half’ indexes i.e., M=min{M: n_(m)=0, |m|∈{M/2, . . . , M}, M=2^(N), N∈

}. Secondly, the distinguished logarithm in (47) is undefined if {circumflex over (γ)}(u)=0. In estimating γ(u) from the data, there is a small but non-zero probability that the estimate will be zero. In this case, the distinguished logarithm in (47) is undefined and the technique fails. As |u| increases, |ϕ_(X)|(u) decreases and may approach |ϕ_(ϵ)|(u). When |ϕ_(X)|(u) and |ϕ_(ϵ)|(u) have similar magnitudes, the probability of |ϕ_(X)+ϕ_(ϵ)| (and hence {circumflex over (γ)}) being close to zero can become significant. The filter H(u) should roll off faster than |ϕ_(X)|(u) approaches |ϕ_(ϵ)|(u) to reduce the impact this may have on the estimate. Ideally H(u) should be zero in regions where noise may result in |{circumflex over (γ)}|(u) being close to zero. Gugushvili has shown [18] that for a rectangular window, the probability of inversion failure approaches zero as the length of the data set increases N→∞.

3.9 Discrete Notation

We digress momentarily to introduce additional notation. Throughout the rest of the paper, bold font will be used to indicate a 2M×1 vector corresponding to a discretely sampled version of the named function, e.g., {circumflex over (ϕ)}_(A) represents a 2M×1 vector whose values are given by the characteristic function {circumflex over (ϕ)}_(A)(u) evaluated at the points u∈{0, 1, . . . , M−1, −M, . . . , −2, −1}. Square bracket notation [k] is used to index a particular element in the vector, e.g., {circumflex over (ϕ)}_(A)[M−1] has the value of {circumflex over (ϕ)}_(A)(M−1). We also use negative indexes for accessing elements of a vector in a manner similar to the python programming language. Negative indexes should be interpreted relative to the length of the vector, i.e., {circumflex over (ϕ)}_(A)[−1] refers to the last element in the vector (which is equivalent to {circumflex over (ϕ)}_(A)[2M−1]).

3.10 Summary of Estimator

The estimation procedure we use may be summarized in the following steps.

-   -   1. Partition sampled time series into intervals using (8).     -   2. Calculate x_(j) value for each interval according to (9).     -   3. Generate histogram n from the x_(j) values.     -   4. Calculate {circumflex over (ϕ)}_(X) using the inverse FFT to         efficiently evaluate (50) at various sample points.     -   5. Calculate ϕ_(Z) and G at the appropriate points.     -   6. Calculate {circumflex over (γ)} via (46) using {circumflex         over (ϕ)}_(X), G and ϕ_(Z).     -   7. Calculate |ϕ_(Asmooth)(u)|, a low-pass filtered version of

${{\frac{1}{\lambda}d{\log\left( \overset{\hat{}}{\gamma} \right)}(u)}}.$

-   -   8. Calculate α_(mod) via (83) and (85).     -   9. Calculate H using α_(mod) and either (86) or (87).     -   10. Calculate {circumflex over (ϕ)}_(A) via (47) using         {circumflex over (γ)} and H. If any element of {circumflex over         (γ)} is zero and the corresponding element of H is non-zero, the         estimation has failed as the distinguished logarithm is         undefined.     -   11. Calculate {circumflex over (ƒ)}_(A) using the FFT of         {circumflex over (ϕ)}_(A) according to

$\begin{matrix} {{{\hat{f}}_{A}\lbrack k\rbrack} = {\frac{1}{2M}{\sum\limits_{m = {- M}}^{M - 1}\;{{{\hat{\phi}}_{A}\lbrack m\rbrack}{e^{- \frac{i\; 2\pi\;{mk}}{2M}}.}}}}} & (88) \end{matrix}$

3.11 Performance Measures

The performance of the estimator is measured using the integrated square of the error (ISE). The ISE measures the global fit of the estimated density.

$\begin{matrix} {{{ISE}\left( {{\hat{f}}_{A},f_{A}} \right)} = {\int_{- \infty}^{\infty}{\left( {{{\hat{f}}_{A}(x)} - {f_{A}(x)}} \right)^{2}{dx}}}} & (89) \end{matrix}$

The discrete ISE measure is given by

$\begin{matrix} {{{ISE}\left( {{\hat{p}}_{A},p_{A}} \right)} = {\sum\limits_{m = {- M}}^{M - 1}\;\left( {{{\hat{p}}_{A}\lbrack m\rbrack} - {p_{A}\lbrack m\rbrack}} \right)^{2}}} & (90) \end{matrix}$

where p_(A) is a 2M×1 vector whose elements contain the probability mass in the region of each histogram bin i.e.,

$\begin{matrix} {{p_{A}\lbrack m\rbrack} = {\int_{m - 0.5}^{m + 0.5}{{f_{A}(x)}{{dx}.}}}} & (91) \end{matrix}$

The vector {circumflex over (p)}_(A) represents the corresponding estimated probability mass vector.

4 NUMERICAL RESULTS OF THE FIRST EMBODIMENT

Experiments were performed using simulated and real data.

4.1 Simulations

The ideal density used by Trigano et al. [11] was used for these simulations. It consists of a mixture of six Gaussian and one gamma distribution to simulate Compton background. The mixture density is given by

$\begin{matrix} {f \propto {{0.5g} + {10{\mathcal{N}\left( {40,1} \right)}} + {10{\mathcal{N}\left( {112,1} \right)}} + {1{\mathcal{N}\left( {50,2} \right)}} + {1{\mathcal{N}\left( {63,1} \right)}} + {2{\mathcal{N}\left( {140,1} \right)}}}} & (92) \end{matrix}$

where

(μ, σ²) is the density of a normal distribution with mean μ and variance σ². The density of the gamma distribution is given by g(x)=(0.5+x/200)e^(−(0.5+x/200)). The density was sampled at 8192 equally spaced integer points to produce the discrete vector p_(A) of probability mass. The FFT was taken to obtain ϕ_(A), a sampled vector of ϕ_(A) values.

-   -   A particular count rate λ was chosen for an experiment,         corresponding to the expected number of events per observation         interval. The expected pile-up density was obtained via (40).         i.e., the discrete vector ϕ_(A) was scaled by λ, exponentiated,         then scaled by e^(−λ) and finally an FFT was applied

$\begin{matrix} {{p_{Y}\lbrack m\rbrack} = {{{{FFT}\left( {e^{- \lambda}e^{{\lambda\phi}_{A}}} \right)}\lbrack m\rbrack}.}} & (93) \end{matrix}$

Equation (93) was convolved with a Gaussian to simulate the effect of noise Z smearing out the observed spectrum

$\begin{matrix} {p_{X} = {p_{Y}\mspace{14mu}\bigstar\mspace{14mu}\frac{1}{\sqrt{2\pi}\sigma}{e^{- \frac{m^{2}}{2\sigma^{2}}}.}}} & (94) \end{matrix}$

This represents the expected density of the observed spectrum, including pile-up and additive noise. Observation histograms were created using random variables that were distributed according to (94). Experiments were parameterized by the pair (N, λ) where N∈{10⁴, 10⁵, 10⁶, 10⁷, 10⁸, 10⁹} and λ∈{1.0, 3.0, 5.0}. For each parameter pair (N, λ), one thousand observed histograms were made. Estimates of the probability mass vector p_(A) were made using (88) with both (86) and (87) used for H(k). A threshold value of α₀=0.95 was used for both window shapes, and β₀=40.0 for the logistic shape. The discrete ISE measure of the error between each estimate {circumflex over (p)}_(A) and the true vector p_(A) were recorded. For comparison with asymptotic bandwidth results, estimates were made using a rectangular window whose bandwidth was selected according to the condition 1.3 specified by Gugushvili in [18] i.e., h_(N)=(ln N)^(−β) where β<½. We emphasize that the β of Gugushvili's filter is not to be confused with the β₀ of (87). The asymptotic bandwidth criterion was implemented by using

$\begin{matrix} {{H\lbrack k\rbrack} = {1_{\{{{k} < \alpha_{0}}\}}\lbrack k\rbrack}} & (95) \\ {{{where}\mspace{14mu}\alpha_{0}} = {\frac{M}{\pi}{\left( {\ln\mspace{14mu} N} \right)^{\beta}.}}} & (96) \end{matrix}$

Three values for Gugushvilli's β were trialed, namely β=½, ⅓, ¼.

-   -   Estimates were also made using a rectangular filter (95) with         fixed bandwidths of various values α/M∈[0.2, 0.4, 0.6, 0.8].         Finally time-series data was created according to (1) with an         idealised rectangular pulse shape and 10⁷ pulses whose energies         were distributed according to (92). The pulse length and count         rate were chosen to give a Poisson rate Δ=1.0. The algorithm         described by Trigano et al. [11] was used to estimate the         underlying amplitude density from a bi-dimensional histogram         containing 32×1024 (duration×energy) bins—this choice of bins         reportedly giving the best accuracy and reasonable execution         times. The performance and processing time of the core algorithm         were recorded for comparison with our proposed algorithm. FIG. 2         plots a typical estimate {circumflex over (p)}_(A) made by the         data-driven logistic shaped filter for an experiment with         parameter pair (N=10⁶, λ=3.0). The true vector p_(A) (thin solid         line), and the observed histogram {circumflex over (p)}_(X)         (lower curve containing some noise) are also plotted. Pile-up         peaks can be clearly seen in the observed histogram. Although         the estimated density suffers from ringing (due to the Gibbs         phenomenon), it otherwise estimates the true density and         corrects the pile-up that was present in the observed histogram.         FIG. 3 plots a typical estimate made at the same operating point         as FIG. 2, but with an estimator having a rectangular filter         where the bandwidth was selected using (96) and β=¼. This         corresponds to the operating region in FIG. 6 where the         performance of the fixed bandwidth filter (β=¼) is approaching         that of the data-driven filters. It is evident that while also         correcting pile-up, the resulting estimate contains more noise.         FIG. 4 shows the distribution densities of ISE measures as a         function of sample count using a rectangular filter and various         fixed bandwidths. Lines were plotted between distribution means         (MISE) to assist visualization. The results for the data-driven         rectangular filter (86) were also plotted, connected with a         thicker curve. This clearly illustrates the weakness of fixed         bandwidth filtering. For any fixed bandwidth, the ISE decreases         as sample count increases, eventually asymptoting as the bias         becomes the dominant source of error. At that point (which is         noise and bandwidth dependent) the ISE remains largely constant         despite increases in sample count. The fixed bandwidth excludes         the use of some estimates {circumflex over (ϕ)}_(A)[k] in the         final calculation, even when they have high         signal-to-noise-ratio (SNR). FIG. 4 also shows the results given         by the rectangular filter with our proposed data-driven         bandwidth selection. This curve lies close to the inflection         point of each fixed bandwidth curve. This indicates the         bandwidth selected for the data-driven rectangular filter is         close to the optimal bandwidth value (for a rectangular filter)         across the range of sample counts. FIG. 5-FIG. 7 show the         distribution densities of the ISE measure as a function of the         total number of estimates N in each histogram at three count         rates λ∈{1.0, 3.0, 5.0}. The MISE curves for the logistic and         rectangular filters are lower than those obtained using the         bandwidth given by (96) for much of the region of application         interest. There are various regions where the non-data-driven         bandwidth (β=¼) gives similar performance to the data-driven         bandwidths, however this is not maintained across the whole         range of sample counts. The logistic filter shape has slightly         better performance than the rectangular filter shape, although         the differences between the two filters appears relatively minor         to the ISE measure. Table 1 compares the results between the         proposed algorithm and the algorithm recently described in [11].         The ISE for both methods were similar at the operating point         under test (λ=1.0, N=10⁷), however our proposed algorithm         requires considerably less computation.

TABLE 1 Comparison With Algorithm Described in [11] Algorithm Avg. ISE Avg. Time (sec) Fast Trigano Algorithm 1.3 × 10⁻⁵ 3.19 32 × 1024 (duration × energy) bins Proposed Algorithm  1 × 10⁻⁵ 0.019

4.2 Real Data

The estimator was applied to real data to assess its usefulness in practical applications. The threshold value ϵ found in (8) was chosen to be approximately one half the standard deviation of the additive noise w(t). This ensured a reasonably high probability of creating intervals, yet ensured errors in the estimation of interval energy were low. A value for the interval length L was chosen approximately four times the ‘length’ of a typical pulse—that is, four times the length of the interval {t: Φ(t)>ϵ}. An energy histogram was obtained from a manganese sample, with a photon flux rate of nominally 10⁵ events per second. A slight negative skew was present in the shape of the main peaks of the observed histogram, suggesting a complicated noise source had influenced the system. This is barely visible in FIG. 8. The noise was modelled as a bimodal Gaussian mixture rather than a single Gaussian peak. A simple least-squares optimization routine was used to fit bimodal Gaussian parameters Z˜α₁

(μ₁, σ₁ ²)+α₂

(μ₂, σ₂ ²) to the noise peak located around bin index zero. A suitable value for λ was chosen manually. The logistic filter with data-driven bandwidth was used to estimate the true density. FIG. 8 shows plots of the observed and estimated probability mass vectors. The main peaks (bins 450˜600) have been enhanced while the pile-up has been attenuated though not fully removed. The first order pile-up peaks have been reduced. The peak-to-pile-up ratio (ratio of the height of main peak to that of first pile-up peak) has increased from around 6 to around 120. These improvements are comparable to other state of the art systems (e.g., [11]). There are several possible reasons the estimator fails to fully resolve pile-up. The accuracy of the estimator depends on correctly modelling the Gaussian noise peak. The bimodal Gaussian mixture modelled the noise peak such that the maximum error was less than 1% of the noise density peak. Given that the residual pile-up peaks in the estimated spectrum are below 1% of the main peak, the sensitivity of the estimator to errors in noise modelling may have contributed to this in some part. A second reason for the unresolved pile-up may be due to the uncertainty in the estimation of the observed spectrum. Several of the residual pile-up peaks are relatively close to the floor of the observed histogram. The residual peaks may simply be a noise induced artefact of the estimator. Finally, the mathematical model may be an overly simple approximation of the observed spectrum. The detection process includes numerous second-order effects that have not been included in the model (e.g., ballistic deficit, supply charge depletion, correlated noise, non-linearities, etc. . . . ). These minor effects may limit the accuracy of the pile-up correction estimator.

5 SUMMARY OF THE FIRST EMBODIMENT

We have taken the estimator proposed by Gugushvili [18] for decompounding under Gaussian noise, and adapted it for correcting pulse pile-up in X-ray spectroscopy. We have proposed a data-driven bandwidth selection mechanism that is easily implemented, and provides significant reduction in ISE/MISE across a broad range of sample counts of interest to spectroscopic applications (10⁴˜10⁹ counts). The data-driven rectangular bandwidth selection is close to optimal (for rectangular filters), and over the range of interest outperforms bandwidth selection based on asymptotic results or fixed bandwidth.

-   -   Although initial results appear promising, further work is         required to improve the performance for practical         implementations. The estimation still contains ‘ringing’         artefacts associated with the Gibbs phenomenon. Additional         filter shape attempts to reduce this, there are other shapes         that are closer to MSE optimal.

6 SECOND EMBODIMENT

This section gives a summary of the spectrum estimator of the 2nd embodiment. The 2nd embodiment solves the problem of the 1st embodiment which requires entire clusters to be approximately encompassed in each interval. In the 2nd embodiment, the entire data series if desired can be used, and the overlap is compensated by introduction of 2 different interval lengths L and L1.

-   -   We need to include a few additional terms not mentioned in the         first embodiment. In particular {circumflex over (ϕ)}_(X1). The         spectrum estimator is based on

$\begin{matrix} {{{\hat{f}}_{A}(x)} = {\frac{1}{2{\pi\lambda}}{\int_{- \infty}^{\infty}{e^{- {iux}}{H(u)}d\;{\log(\gamma)}(u){{du}.}}}}} & (97) \end{matrix}$

The introduction of the filter H(u; α) allows us to address several implementation issues that arise. The estimation procedure we use may be summarized in the following steps.

-   -   1. Partition sampled time-series into fixed length intervals         [T_(j), T_(j)+L), j∈     -   2. Calculate x_(j) value for each interval according to         x_(j)=Σ_(k∈[T) _(j) _(, T) _(j) _(+L))s(k).     -   3. Generate histogram n from the x_(j) values.     -   4. Calculate {circumflex over (ϕ)}_(X) using the inverse FFT of         n.     -   5. Partition the sampled time series into a different set of         intervals with length L₁, and following similar calculations to         obtain {circumflex over (ϕ)}_(X1).     -   6. Calculate ϕ_(Z) and G.     -   7. Calculate {circumflex over (γ)} using {circumflex over         (ϕ)}_(X), G, ϕ_(Z) and {circumflex over (ϕ)}_(X1),

$\begin{matrix} {\hat{\gamma} = \frac{{\hat{\phi}}_{X}}{\phi_{Z}e^{- \lambda}{\hat{\phi}}_{X\; 1}}} & (98) \end{matrix}$

-   -   8. Calculate |ϕ_(Asmooth)(u)|, a low-pass filtered version of

${{\frac{1}{\lambda}d{\log\left( \overset{\hat{}}{\gamma} \right)}(u)}}.$

-   -   9. Calculate α_(mod).     -   10. Calculate H using α_(mod).     -   11. Calculate {circumflex over (ϕ)}_(A) using {circumflex over         (γ)} and H. If any element of {circumflex over (γ)} is zero and         the corresponding element of H is non-zero, the estimation has         failed as the distinguished logarithm is undefined.     -   12. Calculate {circumflex over (ƒ)}_(A) using the FFT of         {circumflex over (ϕ)}_(A) according to

$\begin{matrix} {{{\hat{f}}_{A}\lbrack k\rbrack} = {\frac{1}{2M}{\sum\limits_{m = {- M}}^{M - 1}\;{{{\hat{\phi}}_{A}\lbrack m\rbrack}{e^{- \frac{i\; 2\pi\;{mk}}{2M}}.}}}}} & (99) \end{matrix}$

6.1 Algorithm Details

Partition the detector output stream into a set of non-overlapping intervals of length L i.e., [T_(j), T_(j)+L), T₀∈

_(≥0), T_(j+1)≥T_(j)+L, j∈

_(≥0). Let x_(j) be the sum of the detector output samples in the jth interval i.e.,

$\begin{matrix} {x_{j} = {\sum\limits_{k = T_{j}}^{T_{j} + L - 1}\;{s(k)}}} & (100) \end{matrix}$

Assuming L is greater than a pulse length, the jth interval may contain ‘complete’ pulses as well as pulses which have been truncated by the ends of the interval. It can be shown that x_(j) consists of a superposition of the energy of ‘complete’ pulses which we denote

, the energies of truncated pulses which we denote with

_(1j) and noise z_(j)

-   -   Let the detector output stream be partitioned into a second set         of non-overlapping intervals [T_(1j), T_(1j)+L₁), T_(1,0)∈         _(≥0), T_(1,j+1)≥T_(1,j)+L, j∈         _(≥0) where L₁<L. Let x_(1j) be given by

$\begin{matrix} {x_{1j} = {\sum\limits_{k = T_{1j}}^{T_{1j} + L_{1} - 1}\;{s(k)}}} & (101) \end{matrix}$

If L₁ is chosen to be slightly less than the pulse length, the x_(1j) term will contain no ‘complete’ pulses, but consist of a superposition of only the energies of truncated pulses y_(lj) and noise z_(j). The number of truncated pulses in any interval has a Poisson distribution. We have

$\begin{matrix} {{X_{1} = {Y_{1} + Z_{1}}},} & (102) \\ {\phi_{X_{1}} = {\phi_{Y_{1}}{\phi_{Z_{1}}.}}} & (103) \end{matrix}$

-   -   We can decompose the total energy in the interval [T_(j), T_(j),         +L) into the energy contribution Y₁ from pulses that have been         truncated and the energy contribution Y₀ from pulses that are         fully contained in the interval [T_(j), T_(j)+L) i.e.,

$\begin{matrix} {X = {Y_{0} + Y_{1} + Z_{0} + Z_{1}}} & (104) \end{matrix}$

where Z₀ represents noise in the regions where pulses are fully contained in the interval (a length of L−L₁), and Z₁ represents noise in the regions where pulses are truncated (a length of L₁). Hence,

$\begin{matrix} {\phi_{X} = {\phi_{Y_{0}}\phi_{Y_{1}}\phi_{Z_{0}}{\phi_{Z_{1}}.}}} & (105) \end{matrix}$

-   -   By combining (103) with (105) we have

$\begin{matrix} {\phi_{X} = {\phi_{X_{1}}\phi_{Y_{0}}\phi_{Z_{0}}}} & (106) \end{matrix}$

Rearranging gives

$\begin{matrix} {\phi_{Y_{0}} = \frac{\phi_{X}}{\phi_{X_{1}}\phi_{Z_{0}}}} & (107) \\ {\mspace{76mu}{= {e^{- \lambda_{0}}e^{\lambda_{0}{\phi_{A}{(u)}}}}}} & (108) \end{matrix}$

We can estimate ϕ_(X) ₁ in a similar manner that we estimated ϕ_(X) or some other method, e.g., via the empirical characteristic function or by performing an FFT on the normalized histogram of x_(1j) values.

-   -   When performing the decompounding operation, the Poisson rate λ₀         for the reduced interval length L−L₁ is used to account for the         sub-interval over which the compound Poisson process Y₀ occurs.

6.2 Visualization of Internal Quantities

To aid the reader's understanding, FIG. 9 plots various quantities obtained during the estimation process. The upper blue curve (with a value around 0.3 at bin zero) plots |{circumflex over (ϕ)}_(X)|, the estimated characteristic function of the observed spectrum. A brown curve is used to show the true value of |ϕ_(ƒ)|, which is distinctly visible as the lower curve with periodic nulls in the region [6000, 10000]. The quantity |ϕ_(ϵ)|/(λϕ_(Z)e^(−λ)) is shown in transparent red and appears as ‘noise’ whose average density peaks around bin #8000. The expected value of |ϕ_(ϵ)|/(λϕ_(Z)e^(−λ)) is shown with a black dashed line. This is obtained using (75), the known value of λ, and assuming Gaussian noise with known σ to obtain ϕ_(Z)(k). The quantity |{circumflex over (ϕ)}_(ƒ)| is shown with a transparent blue curve. This is barely visible as it coincides closely with |ϕ_(ƒ)| in the intervals [0, 4000], [12000, 16000], and closely with |ϕ_(ϵ)|/(λϕ_(Z)e^(−λ)) in the interval [5000, 11000]. Note that the colour of |ϕ_(ϵ)|/(λϕ_(Z)e^(−λ)) appears to change from red to purple in the interval [5000, 11000] as both transparent plots overlap. A solid black line shows |{circumflex over (ϕ)}_(fsmooth)|, a low pass filtered version of |{circumflex over (ϕ)}_(ƒ)|. The low pass filtering removes any local oscillations in |{circumflex over (ϕ)}_(ƒ)(k)| due to the peak localities, as described in paragraph on smoothing at the beginning of this section. The term |{circumflex over (ϕ)}_(fsmooth)(k)| serves as an estimate of

[|ϕ_(ƒ)(k)|]. It can be see that |{circumflex over (ϕ)}_(ƒ)| provides a reasonably good estimate of |ϕ_(ƒ)| in the region where |{circumflex over (ϕ)}_(fsmooth)|>>

[|ϕ_(ϵ)|/(λϕ_(Z)e^(−λ))]. As these two quantities approach each other, the quality of the estimate deteriorates until it is eventually dominated by noise. The filter H(k) should include good estimates of |ϕ_(ƒ)| while excluding poor estimates. To find the regions where good estimates of |ϕ_(ƒ)| are obtained we address the question: Given

[|ϕ_(ϵ)/(λϕ_(Z)e^(−λ))], what is the probability that the calculated values of {circumflex over (ϕ)}_(ƒ) in a local region arise largely from noise?

7 COUNT RATE ESTIMATION

The previous estimator assumed λ was known. An estimate of λ can be obtained without prior knowledge as follows.

-   -   1. Using {circumflex over (ϕ)}_(X), {circumflex over (ϕ)}_(X1),         {circumflex over (ϕ)}_(Z), G from the previous section,         calculate

$\begin{matrix} {{\hat{\phi}}_{Y} = \frac{G\;{\hat{\phi}}_{X}}{\phi_{Z}\phi_{X_{1}}}} & (109) \end{matrix}$

-   -   2. Using {circumflex over (ϕ)}_(Y), estimate the count rate.         This can be done a number of ways.     -   3. One way is to use an optimization routine or some other means         to fit a curve to {circumflex over (ϕ)}_(Y). The fitted         parameters can be used to obtain an estimate of the count rate.     -   4. Another way involves estimating the DC offset of Ψ=d         log({circumflex over (ϕ)}_(Y)). This can be done by averaging a         suitable number of points of Ψ. The points obtained by filtering         by H(u) in the previous section are usually suitable, although         less points may also produce an adequate estimate.     -   5. Another way involves using an optimization engine or some         other means to fit a curve to Ψ=d log({circumflex over         (ϕ)}_(Y)). A suitable parameterized curve to fit d         log({circumflex over (ϕ)}_(Y)) is given by

$\begin{matrix} {{f\left( {{u;\lambda},\alpha,\sigma,\mu} \right)} = {{- \lambda} + {\lambda{\sum\limits_{k = 0}^{K - 1}\;{\alpha_{k}{G\left( {\sigma_{k}u} \right)}e^{{- j}\; 2\pi\; u\;\mu_{k}}}}}}} & (110) \\ {{{where}\mspace{14mu}\alpha} = \left( {\alpha_{0},\;\cdots,\alpha_{K - 1}} \right)} & (111) \\ {\sigma = \left( {\sigma_{0},\cdots\;,\sigma_{K - 1}} \right)} & (112) \\ {\mu = \left( {\mu_{0},\cdots\;,\mu_{K - 1}} \right)} & (113) \end{matrix}$

-   -   -   and where K∈             is chosen to allow the curve fit to sufficient accuracy. The             parameter λ provides an estimate of the count rate. The             optimization engine is not required to give equal weighting             to each point in Ψ.

8 DESCRIPTION OF FIGURES

The following figures are to aid understanding of the process. FIG. 1 shows one possible scheme used to partition the detector output. The illustration depicts the sampled detector response to three incident photons. To aid clarity of the figure, the effects of noise have been removed. The output response has been partitioned into several regions of equal length (L). The number of pulses arriving in each region is unknown to the processing system. One pulse has arrived in the first interval. Two pulses have arrived in the second interval. No pulses have arrived in the third interval. The total photon energy arriving within each interval is calculated as the statistic of interest, being the sum of all sample values in each interval. Intervals are not temporally aligned with pulse arrivals. FIG. 2 illustrates the output of the estimation procedure. The true probability density of incident photon energy is plotted as a solid black line. The photons arrival rate is such that three photons on average arriving during any given interval [T_(j), T_(j)+L). The standard deviation of additive noise in the detector output signal s(t) is equal to one histogram bin width. One million intervals were collected. A histogram was made of the total energy in each interval. This is plotted with the blue line. The effects of pile-up is clearly evident, particularly around bins 75, 150 and 225. The red trace plots the estimate of the true incident energy spectrum after the data has been processed by the system. Although some noise appears in the estimate, the effects of pile-up have been removed. The estimate is expected to correctly recover the true incident spectrum on average. This result was obtained using an internal filter whose bandwidth was determined automatically from the data. FIG. 3 illustrates the same quantities as FIG. 2 under the same operating conditions, however in this instance the bandwidth of the internal filter has been determined using asymptotic results from the literature. Although the estimated probability density of incident energies has been recovered, the variance is significantly greater compared to FIG. 2. FIG. 8 illustrates the operation of the system on real data. The blue trace plots the probability density of observed energy values, while the red trace plots the estimated true probability density of incident photon energies. There is no black trace as the true probability density is unknown. In this experiment, X-ray fluorescence of a manganese sample was used as a photon source. The photon arrival rate was around 10⁵ photons per second. The interval length was chosen such that the average time between photons corresponded to the length of two intervals. Sufficient data was collected and partitioned to form 5.9×10⁶ intervals. The standard deviation of the additive noise corresponds to 4.7 histogram bins. The estimation process has clearly reduced the pile-up peaks and enhanced the true peaks. FIG. 9 illustrates various quantities obtained during the simulation of the system described in the 2nd embodiment. It is described in section 5.1 Visualization of Internal Quantities. FIGS. 9-12 relate to the 2nd embodiment. FIG. 10 illustrates the observed and true probability density of input photon energies for the experiment from which FIGS. 9-13 were derived. The black trace plots the true probability density. The red trace plots the density expected to be observed when three photons on average arrive during a given interval length. The blue trace plots the actual observed density. Up to tenth-order pile-up can be seen in the observed density. FIG. 10 included several plots arising from a typical spectroscopic system. The actual incident photon density (‘Ideal Density’) is plotted with a solid dark line. An observed histogram obtained by partitioning the time-series data is shown in dark blue. Distortion of the spectrum caused by pulse pile-up is evident FIG. 13 plots various internal quantities using a logarithmic vertical axis. The dark blue curve that dips in the centre of the plot is |{circumflex over (ϕ)}_(X)|. The green quantity that crosses the plot horizontally is |{circumflex over (ϕ)}_(Y)|. The upper cyan curve that dips in the centre of the plot is |ϕ_(Z)|.1 FIG. 11 illustrates the trajectory of the curve γ in the complex plane. FIG. 12 illustrates internal quantities similar to FIG. 9, however there are some additional signals. The horizontal red trace that is largely noise, and the corresponding black dashed line represent |ϕ_(ϵ)|, the magnitude of the characteristic function of the histogram noise. The transparent green plot that forms the ‘noisy peak’ in the center of the figure is the estimate |ϕ_(ϵ)|. This quantity was plotted in blue in FIG. 9, and was barely visible as it was obscured by |{circumflex over (ϕ)}_(X)|/(λϕ_(Z)e^(−λ)), which is not shown in FIG. 12. The horizontal trace with an average value of −3 is a plot of |ϕ_(Y)|. The cyan trace that begins with a value of zero at bin zero, and dips to a minimum around bin 8000 is |ϕ_(Z)|, the magnitude of the characteristic function of the additive Gaussian noise. FIG. 13 relates to the 2nd broad aspect of calculating the count rate. It illustrates internal quantities used in the calculation of {circumflex over (λ)}. The cyan trace that begins with a value of zero at bin zero, and dips to a minimum around bin 8000 is |ϕ_(Z)|, the magnitude of the characteristic function of the additive Gaussian noise. The dark blue trace that dips to a minimum in the center of the Figure is |{circumflex over (ϕ)}_(X)|, the estimate of the characteristic function of the observed data. The yellow/green horizontal trace with an average value of −3 is the estimate of |ϕ_(Y)|.

REFERENCES

-   [1] G. F. Knoll, Radiation Detection and Measurement, 3rd Edition.     New York: Wiley, 2000. -   [2] P. A. B. Scoullar and R. J. Evans, “Maximum likelihood     estimation techniques for high rate, high throughput digital pulse     processing,” in 2008 IEEE Nuclear Science Symposium Conference     Record, pp. 1668-1672, October 2008. -   [3] M. Haselman, J. Pasko, S. Hauck, T. Lewellen, and R. Miyaoka,     “FPGA-based pulse pile-up correction with energy and timing     recovery,” IEEE Transactions on Nuclear Science, vol. 59, pp.     1823-1830, October 2012. -   [4] T. Petrovic, M. Vencelj, M. Lipoglavsek, R. Novak, and D.     Savran, “Efficient reduction of piled-up events in gamma-ray     spectrometry at high count rates,” IEEE Transactions on Nuclear     Science, vol. 61, pp. 584-589, February 2014. -   [5] B. A. VanDevender, M. P. Dion, J. E. Fast, D. C.     Rodriguez, M. S. Taubman, C. D. Wilen, L. S. Wood, and M. E. Wright,     “High-purity germanium spectroscopy at rates in excess of 10⁶     events/s,” IEEE Transactions on Nuclear Science, vol. 61, pp.     2619-2627, October 2014. -   [6] Y. Sepulcre, T. Trigano, and Y. Ritov, “Sparse regression     algorithm for activity estimation in spectrometry,” IEEE     Transactions on Signal Processing, vol. 61, pp. 4347-4359, September     2013. -   [7] T. Trigano, I. Gildin, and Y. Sepulcre, “Pileup correction     algorithm using an iterated sparse reconstruction method,” IEEE     Signal Processing Letters, vol. 22, pp. 1392-1395, September 2015. -   [8] L. Wielopolski and R. P. Gardner, “Prediction of the     pulse-height spectral distortion caused by the peak pile-up effect,”     Nuclear Instruments and Methods, vol. 133, pp. 303-309, March 1976. -   [9] N. P. Barradas and M. A. Reis, “Accurate calculation of pileup     effects in PIXE spectra from first principles,” X-Ray Spectrometry,     vol. 35, pp. 232-237, July 2006. -   [10] T. Trigano, A. Souloumiac, T. Montagu, F. Roueff, and E.     Moulines, “Statistical pileup correction method for HPGe detectors,”     IEEE Transactions on Signal Processing, vol. 55, pp. 4871-4881,     October 2007. -   [11] T. Trigano, E. Barat, T. Dautremer, and T. Montagu, “Fast     digital filtering of spectrometric data for pile-up correction,”     IEEE Signal Processing Letters, vol. 22, pp. 973-977, July 2015. -   [12] P. Ilhe, E. Moulines, F. Roueff, and A. Souloumiac,     “Nonparametric estimation of mark's distribution of an exponential     shot-noise process,” Electronic Journal of Statistics, vol. 9, no.     2, pp. 3098-3123, 2015. -   [13] P. Ilhe, F. Roueff, E. Moulines, and A. Souloumiac,     “Nonparametric estimation of a shot-noise process,” in 2016 IEEE     Statistical Signal Processing Workshop (SSP), pp. 1-4, June 2016. -   [14] C. McLean, M. Pauley, and J. H. Manton, “Limitations of     decision based pile-up correction algorithms,” in 2018 IEEE     Statistical Signal Processing Workshop (SSP), pp. 693-697, June     2018. -   [15] D. Snyder and M. Miller, Random Point Processes In Time And     Space. New York: Springer-Verlag, 2, revised ed., September 2011. -   [16] B. Buchmann and R. Grübel, “Decompounding: An estimation     problem for Poisson random sums,” Annals of Statistics, pp.     1054-1074, 2003. -   [17] B. van Es, S. Gugushvili, and P. Spreij, “Deconvolution for an     atomic distribution,” Electronic Journal of Statistics, vol. 2, pp.     265-297, 2008. -   [18] S. Gugushvili, Non-Parametric Inference for Partially Observed     Levy Processes. PhD, University of Amsterdam, Thomas Stieltjes     Institute, 2008. -   [19] S. Said, C. Lageman, N. Le Bihan, and J. H. Manton,     “Decompounding on compact Lie groups,” IEEE Transactions on     Information Theory, vol. 56, pp. 2766-2777, June 2010. -   [20] B. van Es, S. Gugushvili, and P. Spreij, “A kernel type     nonparametric density estimator for decompounding,” Bernoulli, vol.     13, pp. 672-694, August 2007. -   [21] J. Yu, “Empirical characteristic function estimation and its     applications,” Econometric Reviews, vol. 23, pp. 93-123, December     2004. 

1. A method of determining a spectrum of energies of individual quanta of radiation received in a radiation detector, the method comprising the steps of: (1) obtaining a time series of digital observations from the radiation detector comprising pulses corresponding to the detection of the individual quanta; (2) computing spectrum sensitive statistics from the detector signal, the spectrum sensitive statistics defining a mapping from a density of amplitudes of the pulses to the spectrum sensitive statistics; and (3) determining the spectrum by estimating the density of amplitudes of the pulses by applying an inversion of the mapping to the spectrum sensitive statistics.
 2. The method of claim 1 further comprising basing the spectrum sensitive statistics on a sum of the digital observations over a plurality of time intervals.
 3. The method of claim 2 further comprising defining the mapping using an approximate compound Poisson process.
 4. The method of claim 3, further comprising augmenting the approximate compound Poisson process by a modelled noise.
 5. The method of claim 4, further comprising expressing the mapping as a relation between characteristic functions of the amplitudes, the spectrum sensitive statistics and the modelled noise.
 6. The method of claim 5, further comprising computing the characteristic functions of the spectrum sensitive statistics by applying an inverse Fourier transform to a histogram of the sum of the digital observations.
 7. The method of claim 5, further comprising computing the characteristic functions of the amplitudes with a low pass filter.
 8. The method of claim 2, further comprising selecting each of the plurality of time intervals to encompass zero or more approximately entire clusters of the pulses, and defining the plurality of time intervals to be nonoverlapping and have a constant length L.
 9. The method of claim 8, further comprising requiring a maximum value of the detector signal at a beginning and end of each time interval.
 10. The method of claim 8, further comprising defining the approximate compound Poisson process as a sum of the amplitudes of the pulses in each time interval.
 11. (canceled)
 12. (canceled)
 13. The method of claim 2, further comprising selecting the plurality of intervals to include: a first set of nonoverlapping time intervals of constant length L without regard to entirety of clusters of the pulses; and a second set of nonoverlapping time intervals of constant length L1 less than L also without regard to entirety of clusters of the pulses; wherein L is at least as long as a duration of the pulses.
 14. The method of claim 13, further comprising selecting L1 to be less than the duration of the pulses.
 15. (canceled)
 16. (canceled)
 17. The method of claim 1, further comprising using a data driven strategy selected to result in a near optimal choice for a kernel parameter which minimizes an integrated square of errors of an estimated probability density function of the energies of the individual quanta of radiation.
 18. A method of estimating count rate of individual quanta of radiation received in a radiation detector, the method comprising the steps of: (1) obtaining a time series of digital observations from the radiation detector comprising pulses corresponding to the detection of the individual quanta; (2) computing spectrum sensitive statistics from the detector signal, based on a sum of the digital observations over a plurality of time intervals, the spectrum sensitive statistics defining a mapping from a density of amplitudes of the pulses to the spectrum sensitive statistics using an approximate compound Poisson process, the plurality of time intervals including: a first set of nonoverlapping time intervals of constant length L selected without regard to entirety of clusters of the pulses; and a second set of nonoverlapping time intervals of constant length L1 less than L also selected without regard to entirety of clusters of the pulses; wherein L is at least as long as a duration of the pulses; (3) determining an estimate of a characteristic function of the approximate compound Poisson process using: ${\hat{\varnothing}}_{y} = \frac{G\hat{\varnothing}x}{\varnothing_{z}{\hat{\varnothing}}_{x}}$ where G is a windowing function, {circumflex over (Ø)}_(X) is an estimate of a characteristic function of a sum of the digital observations over each nonoverlapping time interval in the first set, {circumflex over (Ø)}_(Z) is a characteristic function of a modelled noise process and {circumflex over (Ø)}_(X) ₁ is an estimate of a characteristic function of a sum of the digital observations over each nonoverlapping time interval in the second set; (4) estimating the count rate from the estimate of the characteristic function.
 19. The method of claim 18, further comprising estimating the count rate by using an optimization routine or other means to fit a curve, estimating a DC offset of a logarithm of the estimate of the characteristic function, or fitting a curve to the logarithm of the estimate of the characteristic function. 